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' We study the persistence properties in a simple model of two coupled interfaces characterized 

I by heights hi and h2 respectively, each growing over a d-dimensional substrate. The first interface 

evolves independently of the second and can correspond to any generic growing interface, e.g., of 
the Edwards- Wilkinson or of the Kardar-Parisi-Zhang variety. The evolution of h2, however, is 
coupled to hi via a quenched random velocity field. In the limit d — > 0, our model reduces to the 
Matheron-de Marsily model in two dimensions. For d = 1, our model describes a Rouse polymer 
chain in two dimensions advected by a transverse velocity field. We show analytically that after a 
long waiting time to — > cx3, the stochastic process ^2, at a fixed point in space but as a function of 
time, becomes a fractional Brownian motion with a Hurst exponent, H2 = 1 — /3i/2, where (3i is the 
growth exponent characterizing the first interface. The associated persistence exponent is shown to 

■ be 6's = 1 — H2 = /3i/2. These analytical results are verified by numerical simulations. 
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I. INTRODUCTION 



The survival probability P(t) that a stochastic process X{t) does not cross zero upto time <: is a quantity of long 
. standing interest in probability theory and with many practical applications [1]. The derivative F{t) = —dP/dt is 
^5 ■ called the first-passage probability [2]. This subject has seen a resurgent interest over the last decade in the context 
Q of many body nonequilibrium systems where the stochastic process X{t) denotes a local time dependent field in a 
^ ' spatially extended evolving system. For example, in the case of the Ising model evolving under the Glauber dynamics, 
the relevant stochastic process X{t) is the spin Si{t) at a fixed site j as a function of time t and P{t) then denotes the 
probability that the spin Si does not flip upto time t [3]. In this context, the survival probability P{t) thus measures 
^ , the 'persistence' of a local field to remain in its initial state. In many of these nonequilibrium systems, the persistence 
00 ' has been found to decay as a power law at late times, P{t) ^ t^^ . The exponent 9 is called the persistence exponent 
, and has been a subject of much theoretical, numerical and experimental studies in recent times [4]. The exponent 9 is 
very hard to calculate analytically even in simple systems such as the linear diffusion equation starting from random 
initial conditions [5] . The reason for this difficulty can be traced back to the fact that the spatial interactions in these 
' extended systems makes the local stochastic field X{t) a 'non-Markovian' process in time [4]. 

I Apart from these pure systems, persistence has also been studied in systems with quenched disorder [6,7]. In 
presence of disorder, the exact calculation of persistence P{t) is nontrivial even for a single particle without any 
spatial interaction, though some analytical results have been obtained recently. For example, the asymptotic results 
for the persistence of a particle moving in a random Sinai potential in one dimension have been obtained both in 
the case of a vanishing external field [8-10] and also for nonzero external field [47]. Another solvable example is the 
^ ' persistence of a single particle advected by a layered random velocity field [12-14]. 
O The purpose of this paper is to present analytical results for the persistence P{t) in a system with both spatial 

J-^ interaction and disorder. Our system consists of two growing (rf + l)-dimcnsional interfaces where d refers to the 
dimension of the substrate on which the surfaces grow and 1 refers to the time t. The two interfaces are characterized 
by their heights hi{r,t) and /i2(r,t) respectively. In our model, while the height hi of the first interface evolves 
^ I , independently of /12, the evolution of ft,2 is driven by a velocity that is a random quenched function of hi. This 
' random velocity represents the disorder in the system. The model is detailed in Section II. 
There are two motivations for studying this model: 

(i) We will show later in Section III that our model corresponds to different well known physical systems for different 
values of the spatial dimension d. For example, in the limit d — > 0, our model reduces to the Matheron-de Marsily 
model where one studies the motion of a single Brownian particle in a 2-dimensional plane in presence of a transverse 
random velocity field [15]. On the other hand, for d = 1, we will show that our model describes the evolution of 
a Rouse polymer chain [16] (a Rouse chain consists of a set of beads or monomers connected by harmonic springs) 
moving in a 2-dimensional plane and advected by a transverse random velocity field. While the transport properties in 
this latter model are well understood [17-19], the persistence properties have not been studied before. The analytical 
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results obtained in this paper for the persistence in our model for general d will apply to these models in the limiting 

cases d = and d = 1. 

(ii) The persistence properties of a single interface have been studied extensively both theoretically [20-27] and 
more recently they have been measured experimentally in a system of fluctuating steps on crystal surfaces [28]. In 
the experimental system, there are many step edges, each corresponding to a single interface. If the step edges 
are sufficiently separated from each other, they can be considered as independently growing interfaces. However, in 
general, there will always be an interaction, albeit weak, between two step edges and their motions will be coupled. 
Hence it is important to study persistence in models where the interface heights arc coupled. Motivated by this 
observation, we study here the persistence in a simple model of coupled interfaces. While this model, by no means, is 
a true representative of the actual experimental situation, the advantage is that it represents a minimal model with 
coupling for which one can calculate the persistence properties analytically. 

The paper is organized as follows. In Section II, we define the model precisely and state our main results. In 
Section III we consider the two limiting cases — > and d = 1 where our model reduces to two well known models. 
In Section IV we point out our strategy to compute the persistence exponents analytically. In Section V, we map 
the relevant stochastic process into a fractional Brownian motion (fBm) and then using a well known first-passage 
property of fBm, we calculate the persistence exponents analytically. Section VI considers a generalization to other 
types of interfaces. Section VII contains details of numerical results and finally we conclude with a summary and 
outlook in Section VIII. 

II. THE MODEL AND THE MAIN RESULTS 

We consider two interfaces characterized by their heights hi{r,t) and h2{r,t) growing on a d-dimensional substrate 
according to the following evolution equations 

^=V%+7?i(r,i) (1) 
r)h 

^=V^h2+v{hr{r,t))+V2{r,t), (2) 

where 771 and 772 represent the thermal Gaussian noises that are uncorrelated with each other, each has zero mean 
and their correlators are given by 

{m{ri,ti)mir2,t2)) = [W]-'^/2e-l'-^---l'/4»'j(ti -t2) (3) 
(r?2(ri, ii)r/2(r2, t2)) = [4^a2]-'^/2 e-\-^--\'/^-' d(h - ^2) (4) 

where a denotes the range of the correlator and serves as a short distance cut-off. For d < 2 where the interface 
roughens with time, the short distance cut-off plays no important role and one can safely take the limit o ^ and 
replace the correlators in Eqs. (3) and (4) by delta functions, i.e., [47ra^]^''/^ e^l''i^''^l 6{ri — r2). However, 

for d > 2, the Eqs. (1) and (2) exhibit ultra-violet divergences and one needs to keep a finite a in order that (hi) and 
(/12) remain finite in the t —> 00 limit. 

The first interface hi evolves freely according to Eq. (1) which is precisely the celebrated Edwards- Wilkinson 
equation [29] in {d+ 1) dimensions. On the other hand, the evolution of the second interface /12 in Eq. (2), in addition 
to having the Laplacian and the noise term, is coupled to the height hi of the first interface via the quenched random 
velocity v {hi{r,t)) which is also considered to be a Gaussian with the following moments 



vihi) = (5) 



v{hi)v{h'i) = 6{hi - h[), (6) 



where (. . .) denotes averages over the different realizations of the quenched velocity field v{hi). 
We are interested in the following persistence probabilities at a fixed position r in space 

Pi{t.to) = Prob[/ii(r,t') ^ /ii(r,io) for all t' : to < t' < ta+t] (7) 
P2{t,to) = Prob[/i2(r,i') ^ /i2(r,io) for all t' : to < t' < to + t] . (8) 

The former represents the probability that the height /ii(r,t') of the first interface at a fixed point r in space does 
not return to its value at to within the time interval [^0,^0 + t]- The latter represents the same probability for the 
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second interface. The persistence probability Pi(t, to) for the free Edwards- Wilkinson equation in Eq. (7) has been 
studied both analytically and numerically [20]. It is known that for d < 2, Pi{t,to) has a power law decay for large 
t characterized by a nontrivial persistence exponent and the value of this exponent depends on whether the waiting 
time to = (no return to the initial condition) or to ^ oo (no return to a stationary configuration) [20] 

Pi(i,to = 0) (9) 
Pi(i,to^oo)~t-^^' (10) 

where the subsripts and s in the exponents refer respectively to the probabilities of no return to the initial condition 
(to = 0) and no return to a stationary configuration {Iq oo). The superscript 1 refcrs to the first interface. 
It turns out that the exponent 6^ is hard to determine analytically and is known only numerically. For example, 
^0 = 1-55 ± 0.02 in (1 + l)-dimensions [20]. On the other hand, the exponent 6l was computed analytically [20] by 
mapping the process hi in the to ^ oo limit to a fBm and then using a known first-passage result of fBm. For all 
d <2, one gets [20] 

9l = {2 + d)/4. (11) 

For d = 1 this gives 6] = 3/4, a result that was recently verified experimentally in a system of fluctuating (1 + 1)- 
dimensional steps on Si-Al surfaces [28]. For d> 2, we show that the persistence Pi{t,to — > oo) decays faster than a 
power law for large t, stretched-exponentially for 2 < < 4 and exponentially for d > 4. 

In this paper, we study the persistence probability P2{t,tQ) of the second interface in Eq. (2). As in the case of the 
first interface, we find that for large interval t 

P2{t,to = 0) ^ t-^o (12) 
P2{t,to ^ oo) ^ t-^' (13) 

where the superscript 2 in the exponents refer to the second interface. While we were not able to compute the 
exponent 9q, we calculated the exponent 9^ analytically for all d < 2 

9l = i2-d)/8. (14) 

For d > 2, we argue that the persistence P2{t, to oo) tends to a constant as t ^ oo. 

Furthermore, we were able to generalize this result to the case when the second interface evolves by Eq. (2) but the 
first interface is any generic growing surface, evolving not neccesarily by the Edwards- Wilkinson equation. For example 
the first interface may evolve by the Kardar-Parisi-Zhang (KPZ) equation [30]. In general, this first interface will be 
characterized by a growth exponent f3i and a dynamical exponent zi defined via the scaling form of the second moment 
of the height differences between two points in space, ([fti(ri, n) — /ii(r2, T2)]^) « \t2 — rij^^^/ (jri — r2|^Vk2 — ti|)- 
For example, for the (1-1- l)-dimensional Edwards- Wilkinson equation one has /3i = 1/4 and zi = 2, whereas for the 
(1 + l)-dimensional KPZ equation one has (3i = 1/3 and zi = 3/2 [31]. Our main result is to show that 

O^s = /?i/2. (15) 

In particular, Eq. (15) predicts that in (1-|- l)-dimensions, if the first interface evolves via the KPZ equation, 9^ = 1/6. 
In Section- VII we show that the numerical results are consistent with this theoretical prediction. 

III. LIMITING CASES 

In this Section we show that in the two limiting cases d ^ and d = 1, our model defined by Eqs. (1) and (2) 
reduce respectively to two well studied models. 

A. The Limit d ^ 

In the limit d ^ 0, there is no 'space' variable in the problem. Thus the Laplacian terms on the right hand side 
of Eqs. (1) and (2) drop out and also the noise variables no longer have any r dependence. Interpreting hi = x and 
h2 = y as the coordinates {x,y) of a single particle on a 2-dimensional plane, the Eqs. (1) and (2) reduce to 
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f =,.(«) (16) 
^ = »W<)) + .fe(i). (17) 

These equations represent precisely the MdM model [15] where a single Brownian particle moves in a 2-dimensional 
plane in the presence of a transverse quenched random velocity field v{x). This model was originally introduced to 
study the hydrodynamic dispersion of a tracer particle in porous rocks [15]. While the transport properties in this 
model were well understood [32- 34], the study of persistence properties in this model are relatively recent [12,13]. 
The persistence probabilities defined in Eqs. (7) and (8) reduce, in this context, to the following probabilities 

Pi(t, to) = Prob [x{t') ^ x{ta) for aU t' : to < t' < to + t] (18) 
P2{t, to) = Prob [y{t') ^ y{to) for aU t' : to < t' < to + t] . (19) 

The first probability to) is simply the probability of no return to its initial position of a one dimensional Brownian 
motion and hence Pi{t,to) ~ i"^".' where 9^ ~ Ol = 1/2. The probability P2{t,to) associated with the y coordinate 
is nontrivial. Based on heuristic arguments and numerical simulations, Redner showed that Oq = 1/4 [12]. More 
recently, it was proved analytically that 9o = 9^ = 1/4 hy mapping the stochastic process y{t) to a fBm and then 
using a known first-passage property of the latter [13]. The result 9^ = 1/A is thus a special case of our result in Eq. 
(14) of the interface model in the limit d ^ 0. Incidentally, the exponent Oq = 1/4 happens to be generic for a class 
of transverse velocity fields and occurs even when the tranverse velocity field is a deterministic but an odd function 
of x, i.e., v{x) = —v{—x) [35,36,14]. 



B. The Case d = 1 



Let us consider a Rouse polymer chain embedded in a 2-dimensional plane. The chain consists of beads connected 
by harmonic springs [16]. In addition, the chain is advected by a random layered velocity field as shown in Fig. 1. 
Let [xn{t)-, yn{t)] denote the coordinates of the n-th bead at time t which evolve with time according to the following 
equations of motion 

r(x„+i +x„_i - 2a;„) +r/i(n,t) (20) 

^{Vn+l +yn-l -'2yn) + v{Xn{t)) +r]2{n,t), (21) 

where F denotes the strength of the harmonic interaction between nearest neighbour beads, r]i{n,t) and r]2{n,t) 
represent the Gaussian thermal noises along the x and y directions respectively and are delta correlated. The velocity 
field v{x) is a random quenched function of x taken to be a Gaussian with the following moments 

^x) = (22) 
v{x)v{x') = S{x - x') (23) 

For a finite chain with N beads, the Eqs. (20) and (21) are valid only for the (A'' — 2) interior beads. The two 
boundary beads will have slightly different equations of motion. However, we will only focus here on an infinitely 
large chain (A^ oo) so that the system is translationally invariant along the length of the chain and the boundary 
conditions are irrelevant. 
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FIG. 1. A Rouse chain in a random layered velocity field in 2 dimensions. 



Note that in the absence of the harmonic interaction term, i.e., when F = 0, this model reduces precisely to a single 
particle MdM model discussed in the previous subsection. In presence of the harmonic interaction, the transport 

properties in this model have been studied recently [17 19]. However the persistence properties for F ^ have not 
yet been studied. One can define the following persistence probabilities, 

Pi{t, to) = Prob [xn{t') ^ x„(to) for all t' : to < t' < h + t] (24) 
P2(i, to) = Prob [yn{t') ^ yM for all t' : to < t' < to + t] (25) 

where the former represents the probability that the x coordinate of a bead (say the n-th bead) does not return to 
its position at time to within the time interval [to, to + 1], while the latter represents the same probability for the y 
coordinate of the same bead. For an infinite chain, the system is translationally invariant along the length of the 
chain and hence these persistence probabilities do not depend on the bead label n. 

Since we are interested in the late time properties, one can conveniently replace the harmonic interaction term in 
Eqs. (20) and (21) by a continuous Laplacian operator 

(26) 

^=g^+v{x{s,t))+r]y{s,t), (27) 

where s denotes the distance along the chain and we have rcscaled the time to set the coefficient in front of the 
Laplacian to be unity. Interpreting x = hi and y = h^, the equations (26) and (27) reduce precisely to the d = 1 
version of our model defined in Eqs. (1) and (2). Thus, substituting = 1 in our results for the persistence exponents 
in Eqs. (11) and (14), we obtain the persistence exponents for the polymer problem. In particular, for the probability 
of no return to a stationary configuration [to oo limit in Eqs. (24) and (25)) of the polymer, we have for large 
interval t, 

Pi{t,to ^ oo) ^ t-^/^ (28) 

P2(^,^o^oo)~^-l/^ (29) 

where the exponents 3/4 and 1/8 are obtained by substituting = 1 in Eqs. (11) and (14) respectively. Interestingly, 
for F = 0, i.e., in the single particle MdM model, the corresponding probabilities decay as, Pi(t, to oo) ^ t^^^^ and 
P2it,to oo) ~ t^^/*, as discussed in the previous subsection. Thus, switching on the harmonic interaction strength 
F has opposite effects on Pi and Pj. While the probability Pi{t,to oo) decays faster as ~ t~^/^ (compared to 
~ when F = 0) in presence of interaction F, the probability P2{t,to oo) decays slower as ~ t when F is 

switched on (compared to ~ i^^/^ when F = 0). 



IV. THE STRATEGY TO COMPUTE THE PERSISTENCE EXPONENTS 

The stochastic process hi (r, t) in Eq. (1) for a fixed r, as a function of time t, is Gaussian albeit non-Markovian. The 
non-Markovian property arises due to the Laplacian term that generates interaction between heights at two different 
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space points. The process ft,2(r,t) in Eq. (2) for fixed r is similarly non-Markovian and moreover it is non-Gaussian 
due to the quenched velocity field in Eq. (2). Analytical calculation of the persistence exponent is known to be very 
hard for a non-Markovian process even if the process is Gaussian [4] . For a non-Gaussian process it is even harder in 
general. However, in certain cases, it may be possible to map the relevant stochastic process into a fBm and then one 
can use a known first-passage property of fBm to calculate the persistence exponent analytically. This strategy has 
been successful in the past to compute analytically the persistence exponents in a number of problems even though 
the relevant processes were non-Gaussian and/or non-Markovian [13,20- 22,26,27]. In particular, this strategy was 
used in Ref. [20] to calculate the exponent 9l = (2 -|- (i)/4 in Eq. (11) for the first interface evolving freely with a 
Edwards-Wilkinson dynamics. Here we exploit the same strategy to compute the exponent 9"^ = (2 — d)/8 in Eq. (14) 
for the second interface. 

Before proceeding further, it is then useful at this point to summarize the definition and the known first-passage 
property of a fBm. A stochastic process X{t) (with zero mean E[X{t)] = 0) is called a fBm if its incremental two-time 

correlation function C(ti ,t2) = E {X{ti) — X{t2))'^ is (i) stationary, i.e., depends only on the time difference 1*2 — ti \ 

and (ii) grows asymptotically as a power law 

C(ii,i2)~ |t2-ii|^^, 1*2 -til »1. (30) 

The parameter < if < 1 is called the Hurst exponent that chracterizes the fBm [37] and E[. . .] denotes the 
expectation value over all realizations of the process X(t). For example, an ordinary Brownian motion corresponds 
to a fBm with H = 1/2. The zero crossing properties of a fBm has been studied extensively in the past [38 -40,20]. 
The particular property that is useful for our purpose is the fact that the probability that a fBm does not cross zero 
upto time t decays as a power law at late times, P(t) ~ t~^ with = 1 — H [40,20]. Note that this result 9 = 1 — H 
holds for any zero mean process that satisfies the defining characteristics (i) and (ii) above of a fBm. In particular, it 
holds even if the process is non-Gaussian and/or non-Markovian as long it satisfies (i) and (ii) above. A simple proof 
of this result is given in Ref. [20]. Our strategy in the next Section would be to first show that indeed for d <2, the 
processes /ii(r,i) and /i2(i", *) in Eqs. (1) and (2), for fixed r and in the hmit to oo, do satisfy the characteristics 
(i) and (ii) of a fBm with Hurst exponents Hi = (2 — (i)/4 and H2 ~ {6 + d)/8 respectively and then use the powerful 
result 9 — 1 — H stated above to prove that 6*^ = (2 -I- d)/4 and 9^ = {2 — d)/8. For the free interface hi in Eq. (1) 
this mapping to the fBm was already done in Ref. [20], but we will include it Section V-A for completeness. The 
corresponding calculations for /12 are new and will be detailed in Section V-B. 



V. TWO-TIME CORRELATION FUNCTIONS 



In order to make use of the first-passage property of fBm mentioned in the Section-IV, we need to calculate the 
incremental correlation functions, 



Ci{ti,t2, to) = E [(/ii(r, to + ti) - /ii(r, to + t2)Y 
C2{ti,t2,to) = E f(/i2(r,*o + ii) - h2ir,to+t2)f 



(31) 
(32) 



where E [. . .] denotes an average over both thermal noises as well as the disorder. Note that while for the first interface 
hi, the expectation [. . .] is only over the thermal noise, for the second interface it includes averages over both the 
thermal noise as well as the disorder. In Section V-A, we calculate Ci(*i,*2,*o) which is rather straightforward. The 
computation of C2{ti,t2, to) is more involved and is detailed in Section V-B. 



A. Incremental correlation function for the first interface 



Since the Eq. (1) is linear one can solve it by the standard Fourier transform technique. We define the Fourier 
transform, /ii(k, *) = / d''r /ii(r, *) exp(ik.r). Taking the Fourier transform of Eq. (1) one gets 

hi{h,t) = e"^^' e'''*'fji{k,t')dt', (33) 
Jo 

where we have assumed that the system starts from a flat initial condition hi{r, * = 0) = 0. Using the noise correlator 
in Eq. (3) one can then easily show that 
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The auto-correlation function in real space can then be obtained as 



(34) 



Ai{ti,t2) = {hi(rJ,i)hi{r,t2)) 



I 



-(1 



|f2-tl| 



) + (l-e- 



-/C^(tl+t2) 



(35) 



The right hand side of Eq. (35) is a sum of two integrals each of which is convergent and can be easily done in closed 
form and one gets [20] 



Aiitl,t2) 



ao 



(2-d) 



ih +t2 + a 



2\1- 



{\t2-h\+aY-^ 



(36) 



where ao = (An)-'^/'^. Thus, as t -> oc., Ai{t,t) = {h\{Y,t)) - t^-'^/'^ t^/?! f^j. ^ < 2. One the other hand, for d > 2, 
A\{t,t) ao/{{d— 2)a'^~^) as i ^ 00. Thus while for d > 2 the surface becomes smooth in the stationary state, for 
d <2 the surface roughens with time and the fluctuations grow as a power law in an infinite system 



(/if(r,t)) - t^'^S where /3i = (2 - d)/4. 
The incremental two-time correlation function defined in Eq. (31) can then be writen as 
Ci{ti,t2,to) = {[hi{r,ti+to) - hi{r,t2+to)f) 

= Ai{ti + to,ti+to) + Ai{t2+to,t2+to) - 2Ai{ti + to,t2+to). 



(37) 



(38) 



Using Eq. (36) and taking the limit to — * 00 one finds that for < 2 



Cl{ti,t2,tQ 00) 



{2-d) 

\t2-tl\ 



2(1*2- ill + a 



2\1- 



2a 



2-d 



for 1*2 - hi » 



(39) 



Comparing with the defining property of the fBm in Eq. (30) we see that for d < 2 and in the limit Iq — > 00, i.e., 
when one observes the stochastic process hi{r,t) at a fixed point r in space as a function of time after a waiting 
time to ^ 00, the process is a fBm with a Hurst exponent i?i = (2 — d)/4. This then proves that the persistence 
Pi{t, to 00) ~ t~^' as i ^ 00 where 



el = 1- Hi = {2 + d)/4. 



(40) 
3/4, 



In the limit d ^ 0, Eq. (40) gives 9l = 1/2 (the classical Brownian motion result) and for d = 1, one gets 6l 
a result that was first predicted in Ref. [20] and was subsequently verified experimentally [28]. 

We now turn to the persistence probability Pi{t,to 00) for d > 2. For d > 2, the physics of the process is 
rather different. The fluctutaions of height hi no longer grow at late times, but rather saturates to a constant, 
(/if (r, t)) = Ai{t, t) ao/((d— 2)a''~2) ast^oo. Thus, for d> 2, the cut-off a > is essential. Besides, the relevant 
process hi{r,t), at a fixed r but as a function of time, no longer satisfies the properties of a fBm. Thus one can no 
longer use the first-passage property of fBm to compute the persistence Pi(t,to 00). To make progress, we note 
that for d > 2, it follows from Eq. (36) that in the limit when both ti and t2 are large but their difference 1*2 ^ I = t 
kept fixed, the auto-correlation function becomes stationary, i.e, only a function of * = |*2 — ti\ and decays as a power 
law 



Aliti,t2) 



0,0 



(d-2) ^\t2-ti\+a^f'-'' 



(41) 



Thus, in this hmit, hi is a Gaussian stationary process with an auto-correlator Ai(|*2 — *i| = t) ~ f-i'^/'^-^) , 
Interestingly, exactly the same Gaussian stationary process also represents the total magnetization of a manifold 
in nonequilibrium critical dynamics within the mean field theory [41]. In general, the calculation of persistence of a 
Gaussian stationary process with an algebraically decaying correlator is nontrivial. However, it was pointed out in Ref. 
[41] that one can use a powerful theorem due to Newell and Rosenblatt [42] to obtain useful bounds on the persistence 
property of such a Gaussian process. The Newell-Rosenblatt theorem states that if the stationary auto-correlator of 
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a Gaussian process decays algebraically as A{t) ^ t "for large time difference \t2 — ti\ = t, then the persistence P(t) 
(i.,e., the probability of no zero crossing over a time interval t) has the following asymptotic behaviors, 

P{t)-^ exp[~Kit] for a > 1 
expl-JsTat" Int] < P{t)< exphifgi"] for < a < 1, (42) 

where Ki, and are some constants. In the borderline case a = 1, one has additional logarithmic correction. 
Applying this theorem to our interface problem upon identifying a = (i/2 — 1>0, we conclude that for large t 

Pi{t, to oo)^ e-x.p[-Kit] for d> 4 
exp[-K2t'^/^-^ \nt] < Pi{t, to ^ oo)< exp^ifai'^/^-i] 2 < d< 4. (43) 

Thus the persistence decays exponentially at late times for d > 4 and stretched-exponentially for 2 < d < 4. This is 
to be contrasted with the power law decay Pi{t,to — > oo) for d <2. 



B. Incremental correlation function for the second interface 



In this subsection we show that in the limit to ^ oo and for d < 2, even is a fBm process and one can calculate 
exactly the corresponding Hurst exponent i?2 and hence the persistence exponent 6l = 1 — H2. Taking the Fourier 
transform of Eq. (2) one gets 



/i2(k,t) =6-*='* / dt'e'''*' 
Jo 

Alternatively in real space, one can write 

/i2(r,t) - 



ik.r 



vihiir',t'))d'^r' 



p-(r-r')V4(«-t') 

dt' I „,„ [7^2{r',t')+v{h,{r',t'))]. 



[4n{t - t')]d/^ 

Then the auto-correlation function, averaged over both the thermal noise and the disorder, is given by 



A2{ti,t2) = E[h2{r,ti)h2{r,t2)] = {h2{TM)h2{T,t2)) = HtiM) + HUM), 



(44) 



(45) 



(46) 



where Ii{t\,t2) counts the contribution due to the noise- noise correlator (?72?72) and has the same expression as the 
auto-correlation function of the first interface. 



h{tt,t2) = Ar{tt,t2) 



ao 



(ti +t2 + a^y-^ - i\t2 - ill + ay-^ 



(2-d) 

On the other hand I2{ti,t2) counts the contribution due to the random velocity term and is given by 



dn I dT2 / d'^v[ 



3-(r-r'i)V4(ti-Ti) 

[47r(ii -ri)]'^/2 



d^r: 



g-(r-ri)V4(t2-T2) 
[47r(i2 - T2)]'^/2 



(v.v) 



(47) 



(48) 



where 



{v.v) = {v[hi{r[,Ti)]v[hi{r'2,T2)]) 
= {S{h,ir[^n)-h,{r'2,T2))) 

^ ^gi9[/il(r'i,ri)-?n(ri,r2)]^_ 

, 27r 



(49) 



In Eq. (49) we have first performed the disorder average which gives the expression on the second line. Next we used 
an integral representation of the delta function in the third line. The thermal average over the delta function can be 



done by noting that for a Gaussian process h, {6{h)) = (||)(e*'''*) JZ^ W'f^'^ '^^ ^'"^ 
case h\ being a Gaussian process, (/ii(r'j,Ti) — /ii(r2,T2)) is also Gaussian. This gives 

1 



{v.v) 



V27r((/ii(ri,Ti)-/ii(r^,r2))2) 



l/v/27r(/i2); in our 
(50) 
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Using Eq. (34) one can easily compute the following two-time correlation function 



{hi{r[,Ti)hi{r'^,T2)) 



^2^^^ (/ii(k,ii)/ii(-k,t2))e^'^-("-i-'-^) 



Ti+T2+a 



dT- 



\d/2 ■ 



'|Ti-r2|+a2 (47rT) 

For d <2, one can safely set the cut-off a — > for simplicity and using Eq. (51) one gets 

rn+T2 



(51) 



{[h,{r[,n)-h,{r'„T2)r) 



ao 



(2-d) 



(2ri)i-U(2r2)i-*-(2-ci) / 



\ti—T2\ 



(52) 



Eqs. (52) and (50) thus fully specify I2{ti,t2) in Eq. 48 which can be further simplified by making the change of 
variables: r'^ — r = u and r'2 — r[ = z. The integral over u is a Gaussian that can be easily performed. After a few 
steps of elementary algebra we get 



/2(ii,*2) = V2(2 - d) (47r)('^-2)/4 / d^z 



g-z2/4(ti+t2-Tl-T2) 

[47r {h +t2-Tl- T2)f^ 
ri+T2 



/ dn / dT2 (2ti)^"5 + (2r2)^"^ - (2 - rf) / dr 

Jo Jo J\ti-T2\ 



-I -1/2 



g-zV4r 



(53) 



Note that when ti = t2 = t a simple power counting in Eq. (53) shows that l2{t,t) ~ ^('^+6)74 for large t. On the 
other hand, it follows from Eq. (47) that Ii{t, t) ~ i(2-<^)/2 for large t. Thus l2{t, t) grows faster than Ii{t, t) for large 
t. Thus it follows from Eq. (46) that 



E[hl{Y,t)\ = A2{t,t) t^^\ where /32 = (6 + rf)/8. 



(54) 



In particular for = 1, the result (i2 = 7/8 agrees with that of Ref. [19] derived in the context of a Rouse chain in a 
random tranverse velocity field. 

The incremental correlation function defined in Eq. (32) can then be written in terms of the auto-correlation 
function 



C2{tl,t2,to) = E [(/i2(r,ii +to) - h2{Y,t2+to)) \ ) 

= A2{ti +to,ti+ to) + A2{t2 + to,t2+ to) - 2A2{ti +to,t2+ to) 
= Ci{ti,t2,to) + C^^ , 



(55) 



where in going from the second to the third line in Eq. (55) we have used Eq. (46), the fact that Ci{ti,t2,to) = 
Ii{ti + to,ti + to) + Ii{t2 + to,t2 + to) — 2/1 (ti + to,t2 + to) and the contribution coming from the I2 terms is 
defined as 



Ct = hitl + to, ti + to) + I2{t2 + to, t2 +to)-2 hiti + to, t2 + to). 



(56) 



A close inspection of Eq. (53) reveals that for large ti and t2 the dominant contribution to the integral comes from 
the region where {ti +t2 — tx — T2) 0. To capture the leading behavior of I2{ti,t2) for large ti and <2 it is then 
sufiicient to make the replacement 6'''''/'^^*'+*^-'^'-''^^ /[4:Tr{ti + i2 - n - T2)f^'^ 5{z) in Eq. (53). Thus to leading 
order for large ti and t2 one gets 



/2(ti,i2)~ V2(2-d) (47r)('^-2)/4 



/"*■/' 

Jo Jo 



dT2 



(2ri)i-'^/2 + (2t2)1-''/2 - 2(ri +T2)i-'*/2 + 2|ti - T2\'- 



-d/2 



-1/2 



(57) 



For large to, all the arguments of I2 on the r.h.s of Eq. (56) are large even though ti and ^2 may not be large. Hence 
one can use the asymptotic expression of I2 in Eq. (57). Substituting Eq. (57) in Eq. (56) and rearranging the 
domains of integration one finds for large to 
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« V2(2-rf) (47r)(''-2)/* / rfn / dT2 (2ti)1-<^/2 + (2t2)1-'^/2 - 2(ti + ra)!-'^/^ + 2\n - r^X^-'"'' 

Jtl+to Jti+to '- 



-1/2 



(58) 



Making a further change of variables xi = ti — {ti + to) and X2 = T2 — (^2 + to) one finds that to leading order for 
large to the r.h.s of Eq. (58) becomes independent of to and depends only on \t2 — ti\, 

r\t2-ti\ p\t2-ti\ 

Ci^ « V(2 - d) (47r)('*-2)/^ / da;i / dx2\ Xi - X2 

Ja Jo 

Performing the double integral in Eq. (59) one finally gets 



Substituting the results from Eqs. (39) and (60) on the r.h.s of Eq. (55) one gets for large to 

C2{tut2,to ^ OO) « bl\t2 - tll'-'^/^ + b2\t2 - (61) 

where 61 and 62 are two constants that can be read off Eqs. (39) and (60) respectively. Since the exponent (6 + rf)/4 > 
(1 — d/2) one gets for large 1*2 — ti\, 

C2{tut2,to ^ 00) ~ \t2 - (62) 

Comparing with Eq. (30) one thus finds that for large to and d <2, the process /12 is also a fBm with 

i?2 = (6 + d)/8 (63) 
02 = l-if2 = (2-rf)/8 = /3i/2. (64) 

In the limit d ^ 0, one thus recovers the result 0^ = 1/4 for a single particle MdM model [12,13]. For d = 1, we get 
0'^ = 1 /8, the exact persistence exponent for the Rouse chain in a transverse velocity field. 

The above results are valid for d < 2. For d > 2, one needs to keep the cut-off a finite. Carrying out a similar 
analysis as in the d < 2 case but keeping a finite, one can show that in the limit to 00, the incremental correlation 
function has the following asymptotic behavior for all > 2, 

C2(ii,t2,to^oo)~ |t2-tiP, for 1*2 -til » a^. (65) 

Thus for d > 2, /i2 is a fBm with H2 = I which indicates that 0^ = for all d > 2. This means that the persistence 
P2{t,to — > 00) tends to a constant at large time t for all > 2. Indeed, Eq. (65) indicates that the surface /i2 grows 
ballistically at late times and with a finite probability /12 docs not return to its starting position at to over the time 
interval [*0) ^ + ^o]- In combination with the result in Eq. (64), one thus concludes that the exponent 0^ = {2 — d)/8 
for d < 2 tends to as d — > 2~ and then sticks to 0^ = for all d > 2. Note that this behavior of P2{t,to — > 00) for 
d > 2 is quite opposite to the corresponding d > 2 behavior of the first interface hi for which Pi {t, to — > 00) decays 
faster than a power law at large t as was shown in the previous subsection. 



VI. GENERALIZATION TO OTHER GROWING INTERFACES 



In this Section, we consider a generalization of our model of coupled interfaces. In this generalized version, while 
the second interface height /12 still evolves via Eq. (2), the first interface height hi may correspond to any generically 
growing interface, not necessarily evolving via the Edwards- Wilkinson equation (1). For example, hi may evolve via 
the KPZ equation [30] 

^=^^hi + X{Vhi)'+Viir,t). (66) 
at 

In general, we will consider a generically growing interface hi characterized by the following dynamical scaling of its 
space-time correlation function [31] 
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{[h,ir[,n) - h,{r'„r,)r) « \r, - nf^^f (^J-^) (67) 

where /?i > is the growth exponent, Zi is the dynamical exponent and f{y) is the dynamical scaling function which 
approaches a constant as y ^ 0, /(O) = C and decays for large y. For example, for the (1 + l)-dimensional Edwards- 
Wilkinson equation, /3i = 1/4 and Zi = 2 whereas for the (1 + l)-dimcnsional KPZ equation, f3i = 1/3 and Zi =3/2 
[31]. Note that for a (1 + l)-dimensional KPZ equation, one looses the symmetry h —h. Hence, the persistence 
exponent associated with hi will be different depending on whether the process hi stays above its mean value or 
below its mean value [21,26]. In this paper, wo focus only on the second interface /12 for which there is only one 
persistence exponent 9^. The goal of this section is to show that quite generically 0g = /?i/2. 

We follow a similar route as in Section V-B and start with the calculation of the autocorrelation function of the 
/i2. All the steps from Eq. (44) to Eq. (49) remain unchanged since one doesn't use any information about hi till 
Eq. (49). One uses the specific information about hi for the first time in evaluating the average {v.v) in Eq. (50). In 
Section V-B, the process hi is Gaussian at all times since it evolves via the linear Edwards- Wilkinson equation (1). 
This fact that hi is Gaussian was used explicitly in evaluating the thermal average (exp [i q [/ii(r']^, ri) — hi{r'2, T2)]) in 
Eq. (49) which led to the result in Eq. (50). For a generic non-Gaussian interface hi one can use this step to evaluate 
the thermal average in Eq. (49). To make progress, let us first denote h = [hi{r[,Ti) — /ii(r2,T2)]. Then Eq. (49) 
gives 



where P{h) is the normalized probability distribution of the variable h. For a generic interface, one expects the 
normalized distribution to have the scaling form, P{h) = ^^f^2^ ^ ( \7^^) ' ^^^^^^^^^^'^S ^^^^ scaling form in Eq. 
(68) and rescaling, one finds the following scaling 



(v.v) ^= = =. (69) 

x/((/ii(r;,Ti)-/ii(r^,T2))2) 

Next we substitute the generic dynamical scaling form in Eq. (67) for the correlation function in the denominator of 
Eq. (69) and use the resulting expression of {v.v) on the r.h.s of Eq. (48). The subsequent evaluation of the integral 
-^2(^1,^2) for large ti and t2 followed by the evaluation of the incremental correlation function (^2(^1,^2,^0) can be 
done by following an identical analysis as in Section V-B which we do not repeat here. After performing these steps 
one finds that for large Iq ^ oo 

C2{tl,t2,to ^ 00) ~ 61 1^2 - t2\^-''^^+Ci\t2-ti\''-^\ (70) 

where bi and ci are unimportant constants. Thus, /12 satisfies the defining property in Eq. (30) of a fBm with a Hurst 
exponent given by 

if2 = max[(2-d)/4,l-/3i/2]. (71) 

In particular, for /3i < 1 + d/2 (which seems to be the case for most interfaces), one gets H2 = 1 — pi/2. This then 
leads to the persistence exponent 

el = l-H2= (3i/2. (72) 
For example, for a (1 -|- l)-d KPZ interface hi for which /Ji = 1/3, one gets 

H2 = l- /3i/2 = 5/6; 61 = 1/6. (73) 
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VII. NUMERICAL SIMULATIONS 



In this Section we numerically verify some of the analytical predictions made in the previous sections for the 
coupled interface model, in particular the fact that in the to — > oo limit /12 is generically a fBm with Hurst exponent 
i?2 = 1 - /3i/2, i.e., for large \t2 -ti\, 

~|i2-ti|'-^^ (74) 

where /3i is the growth exponent of the first surface hi . We have checked this prediction numerically in d = 1 for two 
cases: (i) when the first surface evolves via the Edwards- Wilkinson equation so that /?i = 1/4 and Eq. (74) predicts 
C2(ti,t2,U) 00) ~ 1^2 — tiT^'^ for large \t2 — ti\. This case corresponds also to the Rouse chain advected by a 
transverse velocity field as mentioned in Section III-B and (ii) when the first surface evolves via the KPZ equation so 
that Pi = 1/3. In this case, Eq. (74) predicts C2{ti,t2,to 00) ~ \t2 — h}^^^. The results for the simulations in the 
two cases are shown in Figs. (2) and (3) respectively. 

Our simulation techniques are straightforward. For the case (i) above we use the time discretized version of the 
Rouse chain model, i.e., Eqs. (20) and (21). We set tm = mAt and rewrite Eqs. (20) and (21) as 

(U,)] + VAir]i{n,tm) (75) 

yn{tm+l) = yn{tm) + At [yn+l{tm) + Vn-litni) - Vnitm)] 

+ Atv{xn{tm)) + VAir]2{n,t„,). (76) 

For the boundary points n = 1 and n = N,we use free boundary conditions, i.e., we hold a;o = xi, yo = yi, xn = xn+i 
and yN = yw+i for all times tm- We choose At < 0.5 in our simulations so that the stability is guaranteed [20]. The 
variables rii{n,t„i) and 772 (n,m) arc independent Gaussian variables for all n and f„, and each is distributed with zero 
mean and unit variance. We have checked that even if the noise variables have a binary distribution (i.e., +1 and — 1 
each with probability 1/2), the results at long times do not change. Besides, as it turns out from Eq. (70) that the 
thermal noise 772 is actually irrelevant for the long time properties of the /12 process, we have dropped 772 in Eq. (76) 
in most of our simulations. 

We choose the random quenched transverse velocity v{x) in the following way. First we consider a grid along the x 

direction with grid spacing Ax. In fact this grid represents the layered structure of the velocity field. At each point of 
this grid we choose independently v{x) = u{x) / ^ Ax where u{x) is a Guassian random variable with zero mean and 
unit variance. Once a set of {v{x)} is thus chosen, they remain fixed at all times during different thermal histories. 
This set {v{x)} constitutes a particular realization of disorder. Finally one performs the disorder average (. . .) by 
averaging over various realizations of the set {v{x)}. Now, Xn in Eq. (75) at any given time is usually a continuous 
variable and may not correspond to a grid point. In fact, in general, Xn will be between two grid points say a;o and 
xq + Aa;. In such a case, we use, as a convention, v{xo) while simulating the r.h.s of Eq. (76). For a fixed realization 
of the disorder {v{x)}, we average over 30 — 40 thermal histories (generated via rji) and then a final average is done 
over 30 — 40 disorder realizations. 

Using this discretization scheme wc have computed the incremental correlation function C2(ti,t2,to) as defined in 
Eq. (74) beyond some large waiting time to (typically to ~ 30000 steps) for different choices of time step At and the 
grid size Aa;. We have checked that the results do not vary much with the time step or the grid size. We find that, 
as predicted analytically, (72(^1,^2,^0 ~^ 00) depends only on the time difference \t2 — ti\ and for large t2 — ti\, C2 
scales as a power law with an exponent ~ 1.75 which is consistent with the theoretical prediction 7/4. A plot of C2 
as a function of \t2 — t\\ is shown in Fig. (2). This confirms numerically the theoretical prediction that indeed /i2 is 
a fBm with the Hurst exponent H2 = 7/8 and hence also confirms that the persistence exponent 9l = 1/8. 

We also numerically computed the incremental correlation function C2{ti,t2,to 00) for case (ii), i.e., when hi 
evolves by the (1 + l)-dimensional KPZ equation. For the (1 + l)-dimensional KPZ equation in Eq. (66), one needs 
to be careful about the discretization of space and time as has been discussed extensively in the literature [43-47] . In 
this paper we used the discretization scheme proposed by Newman and Swift [45] with periodic boundary condition 

hi{Xn,trn+l) = max [/li (x„+i , t„, ) , (.X„_ i , t,„ ) . /(.i (.X„ , f,„)] + VAir]i{Xn,trn)- (77) 

We show in Fig. (3) the incremental correlation function (72(1*2 — ti \ = t,to ^ oo) versus t. Evidently C2{t) increases 
as a power law with an exponent w 1.695 which is close to the theoretical prediction 5/3 ~ 1.67. The slight discrepancy 
is due to the system size used by us, for which even the scaling exponent f3i for hi actually is not exactly 1/3 but 
close to « 0.305. We have checked that the exponent systematically approaches the expected value with increasing 
N, and here our quoted value is based on the largest N that we could simulate. 



C2{ti,t2,to ^ oo) = lim E 

to— *oo 



{h2{r,to+ti)-h2{r,to+t2)y 
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10000 

FIG. 2. Log-log plot of the incremental correlation function C2{\t2 — ti \ — t) versus t for the Rouse model. The chain length 

is = 1024, while At = 0.1 and Ax — 0.5; number of disorder histories = 30 and thermal histories = 30. The numerical data 
(shown by + signs) is compared to the theoretical prediction of a power law with power 7/4 as shown by the straight line. 




10000 

t 

FIG. 3. Log- log plot of the incremental correlation function C2(|t2 — ti\ = t) versus t when hi evolves via the 
(1 + l)-dimensional KPZ equation. The chain length is iV = 4096, while At = 0.1 and Aa; = 0.5. The theoretical pre- 
diction is the straight line that corresponds to a power law in t with a power 5/3 ~ 1.67, while a straight line fit to the 
numerical data (shown by -|- signs) gives an exponent w 1.695. This slight discrepancy is commented upon in the text. 



VIII. CONCLUSION 



In this paper we have studied the time-dependent properties in a simple model of coupled interfaces characretized 
by heights hi and /12 respectively growing over a rf-dimensional substrate. The evolution of the first interface hi is 
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not affected by /i2- In fact, hi can be any generic growing interface characterized by a growth exponent For 
example, hi may be evolving either via the Edwards- Wilkinson equation or via the Kardar-Parisi-Zhang equation. 
The evolution of the second interface /i2 however is coupled to that of hi by a transverse quenched random velocity 
field, in addition to having the usual Edwards- Wilkinson dynamics. In the limiting case d ^ 0, our model reduces to 
the Matheron-de Marsily model where one studies the motion of a Brownian particle in a 2-d plane in presence of a 
transverse velocity field. In the limit d= 1, our model describes the equations of motion of a Rouse polymer chain in 
two dimensions in presence of a transverse velocity field. 

We have obtained analytical results for the persistence properties in this model. The main result of this paper is to 
show analytically that after a long waiting time {to — » oo), the process /12, at a fixed point in space but as a function 
of time, becomes a fractional Brownian motion with a Hurst exponent H2 = 1 — /3i/2. By using a known first-passage 
property of fractional Brownian motion we have then shown that after a long waiting time to 00, the persistence 
probability P2{t) that the process /i2 at a fixed point in space does not come back to its value at to over the time 
interval [^0,^0 +t] decays as a power law for large t, P2{t) ^ t~^= where the exponent 9l ~ 1 — H2 = /3i/2. For = 1, 
these analytical predictions have been verified numerically in two cases: when hi evolves via the Edwards- Wilkinson 
equation and when hi evolves via the Kardar-Parisi-Zhang equation. 

The mapping of a relevant stochastic process in some limits to a fractional Brownian motion is a rather powerful 
technique for studying the first-passage properties for non-Gaussian and/or non-Markovian processes. The mapping 
does not work always, but if it works one can use the known first-passage property of the fractional Brownian motion. 
This technique has been used successfully in a number of contexts previously [13,20 22,26,27]. We have demonstrated 
that the same technique also works in another class of coupled interface models discussed in this paper. It would 
be interesting to find other such cases where one can apply the same technique successfully. Finally, it would be 
interesting to study the persistence properties in more realistic models of coupled interfaces that are closer to the 
experimental situation of fiuctuating steps [28] . 
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